# Library -----------------------------------------------------------------

library(mizer)
## Warning: package 'mizer' was built under R version 4.4.3
library(mizerExperimental)
## 
## Attaching package: 'mizerExperimental'
## The following object is masked from 'package:mizer':
## 
##     validSim
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.4.3
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'tibble' was built under R version 4.4.3
## Warning: package 'tidyr' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'purrr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
## Warning: package 'stringr' was built under R version 4.4.3
## Warning: package 'forcats' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Setting parameters ------------------------------------------------------

# Set parameters as a single species model
params_initial <- newSingleSpeciesParams()

# Parameters with a reduced resource level
params_rr <- setResource(params_initial,
                      resource_dynamics = "resource_semichemostat")


# Parameters with initial biomass doubled and initial resources halved
params_double <- params_rr
params_double <- setResource(params_double,
                             resource_level = 0.1)
initialN(params_double) <- 2*initialN(params_double)
initialNResource(params_double) <- initialNResource(params_double)/2

# Fishing gear ------------------------------------------------------------

# Set up fishing gear using sigmoid_weight() selectivity function

gear_params(params_double) <- data.frame(
  gear = "gear",
  species = "Target species",
  catchability = 0.3,
  sel_func = "sigmoid_weight",
  sigmoidal_weight = 15,
  sigmoidal_sigma = 5)

params_double <- setFishing(params_double, gear_params = gear_params)

gear_params(params_double)
##                      gear        species catchability       sel_func
## Target species, gear gear Target species          0.3 sigmoid_weight
##                      sigmoidal_weight sigmoidal_sigma
## Target species, gear               15               5
# Simulation -------------------------------------------------

# Simulate biomass density when initial biomass is doubled and initial resources
# are reduced
sim_double <- project(params_double, t_max = 15, effort = 1)
animateSpectra(sim_double, total = TRUE, power = 2, 
               ylim = c(1e-8, NA), wlim = c(1e-3, NA))
# Shows biomass level oscillating --> predator prey relationship

# Extract yield over time dependent on the fishing gear
getYieldGear(sim_double)
## , , sp = Target species
## 
##     gear
## time         gear
##   0  0.0012712274
##   1  0.0004614314
##   2  0.0002662361
##   3  0.0004220103
##   4  0.0008738575
##   5  0.0007496112
##   6  0.0004361817
##   7  0.0003875088
##   8  0.0005553817
##   9  0.0006736737
##   10 0.0005680980
##   11 0.0004720735
##   12 0.0005009610
##   13 0.0005914342
##   14 0.0006062247
##   15 0.0005487976
plotYieldGear(sim_double)

# Figures -----------------------------------------------------------------


# Figure showing flux over increasing biomass density

# Extract weight and count number
N <- finalN(sim_double)["Target species", , drop = TRUE]
w <- w(params_double)

# Calculate individual growth rate
E_growth <- getEGrowth(params_double)["Target species", , drop = TRUE]
gr <- w * E_growth
flux <- gr * N


flux_data <- data.frame(
  Weight = w,
  Flux = flux
)

# Plot figure
flux_plot <- ggplot(flux_data, 
       aes(x = Weight, y = Flux)) +
  geom_smooth() +
  scale_x_continuous(expand = c(0, 0), limits = c(0,105)) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA)) +
  labs(
    x = paste0("Weight (g)"),
    y = paste0("Flux (g/year)"),
    title = "Flux over increasing weight of fish") +
  theme_classic()

flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_smooth()`).

ggsave("figures/flux_plot.png",
       plot = flux_plot,
       device = "png",
       width = 8,
       height = 6,
       units = "in",
       dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_smooth()`).
# Plot flux on a log-log axis
flux_plot <- ggplot(flux_data, 
                    aes(x = Weight, y = Flux)) +
  geom_smooth() +
  scale_x_log10() + 
  scale_y_log10() +
  labs(
    x = paste0("Weight (g)"),
    y = paste0("Flux (g/year)"),
    title = "Flux over increasing weight of fish") +
  theme_classic()

flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggsave("figures/flux_plot_with_fishing.png",
       plot = flux_plot,
       device = "png",
       width = 8,
       height = 6,
       units = "in",
       dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
# Changing parameters -----------------------------------------------------
# Make local reductions in resource replenishment rate and resource capacity

# Calculate flux
N <- finalN(sim_double)["Target species", , drop = TRUE]
w <- w(params_double)

E_growth <- getEGrowth(params_double)["Target species", , drop = TRUE]
gr <- w * E_growth
flux <- gr * N


# Plot flux before local reductions (without fishing gear)
initial_flux_data <- data.frame(Weight = w, 
                             Flux = flux)
initial_flux_plot <- ggplot(initial_flux_data,
                         aes(x = Weight, y = Flux)) +
  geom_smooth() +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    x = paste0("Weight (g)"),
    y = paste0("Flux (g/year)"),
    title = "Flux over increasing weight of fish") +
  theme_classic()
initial_flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# Reduce rr
rr <- resource_rate(params_double)
w_full <- w_full(params_double)
w_full[215:240]
##  [1] 0.04466836 0.05011872 0.05623413 0.06309573 0.07079458 0.07943282
##  [7] 0.08912509 0.10000000 0.11220185 0.12589254 0.14125375 0.15848932
## [13] 0.17782794 0.19952623 0.22387211 0.25118864 0.28183829 0.31622777
## [19] 0.35481339 0.39810717 0.44668359 0.50118723 0.56234133 0.63095734
## [25] 0.70794578 0.79432823
rr[215:240] <- rr[215:240] / 100000000


# Plot resource rate against weight to check reduction
rr_data <- data.frame(
  Weight = params_double@w_full,
  rr = rr
)

rr_plot <- ggplot(rr_data, 
                  aes(x = Weight, y = rr)) +
  geom_smooth() +
  scale_x_log10() + 
  scale_y_log10() +
  theme_classic()
rr_plot
## Warning in scale_y_log10(): log-10 transformation introduced infinite values.
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_smooth()`).

# Reduce resource capacity
rc <- resource_capacity(params_double)
w_full[215:240]
##  [1] 0.04466836 0.05011872 0.05623413 0.06309573 0.07079458 0.07943282
##  [7] 0.08912509 0.10000000 0.11220185 0.12589254 0.14125375 0.15848932
## [13] 0.17782794 0.19952623 0.22387211 0.25118864 0.28183829 0.31622777
## [19] 0.35481339 0.39810717 0.44668359 0.50118723 0.56234133 0.63095734
## [25] 0.70794578 0.79432823
rc[215:240] <- rc[215:240] / 10000

# Plot resource capacity over weights to see reduction
rc_data <- data.frame(
  Weight = params_double@w_full,
  rc = rc
)

rc_plot <- ggplot(rc_data, 
                  aes(x = Weight, y = rc)) +
  geom_smooth() +
  scale_x_log10() + 
  scale_y_log10() +
  theme_classic()
rc_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

params_reduced <- setResource(params_double,
                              resource_capacity = rc,
                              resource_rate = rr,
                              balance =  FALSE)

gear_params(params_reduced)
##                      gear        species catchability       sel_func
## Target species, gear gear Target species          0.3 sigmoid_weight
##                      sigmoidal_weight sigmoidal_sigma
## Target species, gear               15               5
# Narrow predation kernel
pred_kernel <- getPredKernel(params_reduced)

pred_kernel_reduced <- pred_kernel[, , 215, drop = FALSE]

ggplot(melt(pred_kernel_reduced)) +
  geom_line(aes(x = w_pred, y = value)) +
  scale_x_log10()

select(species_params(params_reduced), beta, sigma)
##                beta sigma
## Target species  100   1.3
given_species_params(params_reduced)$sigma <- 0.2
given_species_params(params_reduced)$beta <- 100

getPredKernel(params_reduced)[, , 180, drop = FALSE] %>% 
  melt() %>% 
  ggplot() +
  geom_line(aes(x = w_pred, y = value)) +
  scale_x_log10()

# lower maximum intake rate
species_params(params_reduced)$h
## [1] 30
given_species_params(params_reduced)$h <- 30


# Calculate new flux
sim_reduced <- project(params_reduced, t_max = 15, effort = 1)

N_reduced <- finalN(sim_reduced)["Target species", , drop = TRUE]
w <- w(params_reduced)

E_growth_reduced <- getEGrowth(params_reduced)["Target species", , drop = TRUE]
grr <- w * E_growth_reduced
flux_reduced <- grr * N_reduced


reduced_flux_data <- data.frame(Weight = w, 
                             Flux = flux_reduced)

reduced_flux_plot <- ggplot(reduced_flux_data,
                         aes(x = Weight, y = Flux)) +
  geom_smooth() +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    x = paste0("Weight (g)"),
    y = paste0("Flux (g/year)"),
    title = "Flux with locally reduced resource replenishment rate and resource capacity") +
  theme_classic()

initial_flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

reduced_flux_plot
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# Save plots
ggsave("figures/flux_plot_initial.png",
       plot = initial_flux_plot,
       device = "png",
       width = 8,
       height = 6,
       units = "in",
       dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggsave("figures/flux_plot_reduced.png",
       plot = reduced_flux_plot,
       device = "png",
       width = 8,
       height = 6,
       units = "in",
       dpi = 300)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'